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We study in this paper the Blume-Emery-Griffiths model in a thin film of stacked triangular 
lattices. The model is described by three parameters: bilinear exchange interaction between spins 
J, quadratic exchange interaction K and single-ion anisotropy D. The spin Si at the lattice site 
i takes three values (—1,0, +1). This model can describe the mixing phase of He-4 {Si = +1, —1) 
and He-3 {Si = 0) at low temperatures. Using Monte Carlo simulations, we show that there exists 
a critical value of D below (above) which the transition is of second-(first-)order. In general, the 
temperature dependence of the concentrations of He-3 is different from layer by layer. At a finite 
temperature in the superfluid phase, the film surface shows a deficit of He-4 with respect to interior 
layers. However, effects of surface interaction parameters can reverse this situation. Effects of the 
film thickness on physical properties will be also shown as functions of temperature. 


PACS numbers: 05.70.Fh, 75.10.-b, 75.70.-i 

I. INTRODUCTION 

The physics of thin films has seen a spectacular de¬ 
velopment during the last 30 years. This is due, on the 
on hand, to numerous electronic applications using thin 
films [l|-|3|, superlattices and multilayers Q, and on the 
other hand, to the lack of theoretical understanding of 
surface properties which are very different from the bulk 
ones. Since it is rather easy to change the conditions at 
the surface of a films, by coating or by adsorption of other 
species, for instance, surface physics offers a lot of oppor¬ 
tunities to discover new microscopic phenomena leading 
to potential electronic applications. One has seen in re¬ 
cent years applications using phenomena such as giant 
magnetoresistance [El, [6|, spin transfer torques Q, spin 
valves, etc.jil. 

Theoretically, surface effects in thin films such as sur¬ 
face phonons, surface magnons, surface plasmons have 
been widely studied. We will focus in this paper on the 
magnetic properties of thin films. In surface magnetism, 
much has been understood, in particular on the existence 
of surface-localized spin-waves and their effects on phys¬ 
ical behaviors of thin films at finite temperatures such 
as the reduction of the critical temperature and the low 
surface magnetization ikH- In most of previous stud¬ 
ies, the spin models such as Ising and Heisenberg models 
have been widely investigated. 

In this paper, we use the Blume-Emery-Griffiths 
(BEG) model to study physical behaviors of thin films. 
The spin in this model has three states +1,-1 and 0. A 
site with a value 0 represents a vacant site. The system 
is considered as a dilute magnetic systems in which the 
number of vacant sites varies as a function of tempera¬ 
ture (T). This model can also describe the mixing phase 
of superfluid He-4 {Si = +1, — I) and normal fluid He-3 
{Si = 0) at low temperatures [l2, [i3- Other models ex¬ 


tended from the original BEG model have been recently 
introduced to study the effects of vacancies and of the 
continuous degrees of freedom in the mixtures He-3 and 

He-4 [3El. 

The present work has been motivated by the desire, 
on the one hand, to know if results of bulk BEG model 
remain valid in films, and on the other hand to see if 
the transition criticality can be altered when we reduce 
the film thickness as we have seen in Ref. |T^ and ES 
(i) there is a cross-over from 3-dimensional (3D) criti¬ 
cality to 2-dimensional (2D) universality with decreasing 
thickness for a second-order transition [TgI, and (ii) the 
3D first-order transition becomes a second-order transi¬ 
tion at very small thickness E3- We note that the BEG 
model has been studied by a number of authors in thin 
films of simple cubic lattice structure [Tsl-I^ but the mo¬ 
tivation was not the same as the present paper. In Refs. 
El and El the authors have used the mean-field approx¬ 
imation to study the phase diagram of a double-layer 
and 5-layer films using negative values of biquadratic K 
term in Eq. (HD below. They found a very rich phase dia¬ 
gram with a tricritical point and a staggered quadrupolar 
phase. The authors of Ref. H studied the same model 
of 5-layer film with negative K by mean-field theory. To 
our knowledge, the mean-field theory cannot be used to 
determine precise phase diagrams and the criticality in 
general in 2D and 3D systems. In Ref. [ 2 TI, the authors 
have investigated a system of mixed spins S = 3/2 and 
2 using the BEG model. They determined the phase di¬ 
agram in various parts of the phase space. Very simple 
Monte Carlo (MC) simulations have been carried out. 
In particular, too few data were obtained in the critical 
phase transition region to be useful for the determination 
of the transition temperature. 

Section II is devoted to the description of the model 
and the simulation method. Results will be shown and 
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discussed in Sect. III. Concluding remarks are given in 
Sect. IV. 

II. MODEL AND METHOD 

A. Model 


B. Method 

We use MC simulation [22| to calculate properties of 
the system at finite temperatures for the size of LxLxLz 
where Lz is the film thickness. Periodic boundary condi¬ 
tions are used in the xy plane. The standard MC method 
is used to study the phase transition. The averaged en¬ 
ergy and the specific heat per spin are defined by 


The Blume-Emery-Griffiths (BEG) model consists of a 
system with three states per spin. The model is described 
by the Hamiltonian 



( 2 ) 


h = -jY. ^ E ^ E (1) 

<i,j> <ij> i 


where the spin variable takes the value Si = —1, 0,1 and 
^<ij> denotes a summation over all nearest-neighbors 
(NN). The model presents in addition to the bilinear spin 
interaction J a biquadratic interaction K and a single¬ 
ion crystal field D. This model can describe the mixing 
phase of He-4 (5,; = +1,-1) and He-3 {Si = 0) at low 
temperatures [l2. fl^. 

The present paper studies this model in a film com¬ 
posed of Lz infinite xy triangular lattices stacked in the 2 ; 
direction. We have chosen the stacked triangular lattices 
to have a larger coordination number of neighbors. Since 
we worked with very thin films of small quantities of mat¬ 
ter, such a large coordination number reduces numerical 
errors on statistical fluctuations. Each site is occupied 
by a spin of values ±1,0. In the Helium language, these 
spins are atoms He-3 and He-4. Since we are working 
at a given finite temperature T (canonical method) we 
leave the system to determine the concentration of He-3 
and He-4 at equilibrium at each given T. The mixing of 
He-3 and He-4 is at random corresponding to the max¬ 
imum entropy. As will be seen below, when there is a 
surface the He-3 concentration is more important near 
the surface than in the interior. 

In the bulk, if D = 0 then the model is Ising-like. 
Nonzero values of chemical potential D favor the prolif¬ 
eration of zero spins in the system and lower the transi¬ 
tion temperature between the superfluid and the normal 
fluid phase. By increasing D the system can support su¬ 
perfluid ordering but with a mixture of the normal liquid 
(He-3) and the superfluid one (He-4). The interplay be¬ 
tween the superfluid-like ordering bilinear term and the 
phase breaking D term generates an exotic phase dia¬ 
gram that consists of a line of continuous phase transi¬ 
tions at low D values and high temperatures [IMl. At 
high D and low T values, the transition becomes discon¬ 
tinuous. 

We will show below that the bulk feature of the phase 
diagram in the space (T, D) is found for thin films though 
there is a variation of the critical value of D above which 
the transition is of first order (see next section). We will 
also show that the phase diagram depends on the surface 
parameters. 


Cv 


= N 


{E^) - {E)^ 

ksT^ 


( 3 ) 


where (...) indicates the thermal average. We define the 
order parameter Q for the g-state Potts model by 

Q = [gmax(( 5 i,( 52 ,Q 3 ) - f]/{q - 1 ) ( 4 ) 

where Qn is the spatial average defined by 


1 ^ 


( 5 ) 


i=l 


n = —1,±1,0 indicates the value of spin Si at site i, 
S{Si^n) the Kronecker symbol, and N the total number 
of sites. Our choice of the 3-state Potts order parameter 
is motivated by the fact that the values of the parameter 
-1, ±1 and 0 make it impossible to use the usual definition 
of the magnetization where one sums all values of spins. 
It is however to be used with complementary quantities 
such as the average number of each of the values 1,-1 and 
0 because the Potts order parameter can tell us if there 
is an ordering but it does not give the information on 
the value of spin which makes the system ordered. The 
susceptibility per spin is defined by 


X 


ksT 


(6) 


In general, we discard about 10^ MG steps per spin 
to equilibrate the system at temperature T before av¬ 
eraging physical quantities over the next 10^ MG steps. 
We shall show below an example of the time evolution of 
the energy and the order parameter. Eor histograms, we 
recorded in general 10^ MG steps per spin. The lattice 
sizes used in our simulations are L = 20, 30,..., 120, 300 
and Lz = 4, 8,12,16. 


III. RESULTS 

Let us take J = 1 and iX = 1, namely ferromagnetic 
interactions between NN. Before showing the results for 
thin films, let us show first the results for the bulk proper¬ 
ties of of the BEG model applied to the stacked triangu¬ 
lar lattices. These results by symmetry argument do not 
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FIG. 1: (Color online) Bulk properties: Energy (top) and 
magnetization (middle) versus T around the tricritical D of 
the bulk case. From right to left: D — 7.3, 7.4, 7.5 and 7.6. 
The cross-over from the second order to the hrst order occurs 
at Dc — 7.5. The transition temperature versus D is shown 
in the bottom hgure: the arrow indicates the bulk tricritical 
point. Simulations have been done for a crystal of 20^ sites 
with periodic boundary conditions in all directions. 


bring new physics with respect to the case of simple cu¬ 
bic lattice [13 • However, these results provide elements 
for comparison with the film case which will be shown 
in details below. We show in Fig. [T]the bulk case: the 
top figure shows the energy E versus temperature T for 
several values of D in the tricritical region. As seen, the 
transition is continuous for D < Dc — 7.5 and discon¬ 
tinuous for D > Dc. The magnetization (middle figure) 
shows also a discontinuity for D > Dc. The bottom fig¬ 
ure shows the transition curve in the space {D^T) where 
the tricritical point is indicated by an arrow. Technical 
details on the determination of Dc are similar to those 
used in the case of thin films. They will be given below. 

For a given film thickness, we study in the same man¬ 
ner the behavior of the BEG model for different values 


FIG. 2: (Golor online) Energy E (top) and magnetization M 
(middle) versus T for D = 6, Lz =4 and L = 120. The 
size effect on E in the transition region is zoomed for L = 36 
(red) and 300 (blue). The size effect on M for those sizes is 
not clearly distinguished. 


of D by calculating the energy, specific heat, the order 
parameter, the layer magnetization and the energy his¬ 
togram. 


A. Order of the phase transition 

In Fig. [2] we show the energy E and magnetization M 
versus T in the case where D = 6, = 4 with L = 120. 

Note that the lateral size effect is slightly seen in E as 
zoomed in the bottom figure, but it is not distinguishable 
in M. The curves E and M present a second-order phase 
transition at Tc — 3.82. 

With increasing D^ the system undergoes a first-order 
transition. We show in Fig. [3] the case of D = 7.3 where 
one observes a discontinuity at the transition tempera¬ 
ture Tc ^ 2.694. 
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FIG. 3: (Color online) Energy (top) and magnetization (bot¬ 
tom) versus T in the first-order region of D\ D = 7.3, 
L = 120, L, = 4. 



FIG. 4: (Color online) Energy histograms at Tc = 3.820 (top) 
and Tc = 2.694 (bottom) for T) = 6 (second-order transition) 
and 7.3 (first-order transition), respectively. 


Using the histogram technique |23l-l25j , we explored the 
transition region to search for the nature of the transi¬ 
tion. For T) = 6, we obtain only a one-peak structure at 
the critical temperature (see Fig. 01 top). The energy 
histogram taken at Tc in the case D = 7.3 exhibits a 
double-peak structure as shown in Fig. 0] (bottom), con¬ 
firming thus the first-order character of the transition 


At a first-order transition the ordered and disordered 
phases coexist. In most cases, the system has mixed 
domains of two phases at the same time, the energy of 
the system is thus the average of the energies of the two 
phases {Ei + £’ 2 )/ 2 . It is however possible that at the 
transition the system goes back and forth between the 
two phases during the time evolution. This is what we 
observe here: we show in Fig. [5] how the energy and the 
magnetization evolve during the equilibrating time of 10^ 
MC steps/spin. There are several remarks: 

(i) In a general manner, in MC simulations a trick to 
use to check the equilibrium time is to do two simula¬ 
tions one with a random initial spin configuration and 
the other one with the ground-state configuration. We 
monitor various physical quantities with time evolution. 
The equilibrium is attained when two initial spin configu¬ 
rations give the same results. We see in Fig. [5] that only 
after a few thousands of MC steps that the two initial 
configurations give statistically the same results. 

(ii) The evolutions of E and M show bimodal distri¬ 
butions over periods of 10^ MC steps. The time of 10^ 
MC steps for equilibrating and 10^ MC steps for averag¬ 


ing is thus sufficient as said above. 

We have calculated the transition temperature with 
varying D from 0 to 7.5. The maximum value of D for 
a 4-layer film is 7.5 above which there is no transition at 
all. This value depends on the film thickness. It comes 
from the fact that the maximum of D should cancel the 
energy from J and K term. For example, with Lz = 4, 
the energy of J and K terms is 

El = — (7J + 1K)2 (2 surfaces) 

— (8J + 8iF)2 (2 interior layers) = —30 

where J = K = 1. The energy from D is E 2 = -\-2D 
(2 surface atoms)+2T) (two interior atoms) = AD. The 
maximum of D is determined by setting Ei T E 2 = 0, 
from which D = 30/4 = 7.5. The same calculation can 
be done for another thickness, yielding another value of 
maximal D. 

To determine the critical value of namely Dc^ where 
the transition changes from second to first order, we fol¬ 
low the variation of the energy gap AT defined by the 
energy separation of the two peaks in the energy distri¬ 
bution. This gap is zero when the transition is of second 
order because the energy distribution is continuous. Us¬ 
ing the histogram method with various values of T, we 
show in Fig. [6] (top) the variation of AT versus D. As 
seen, AT is not zero for D g]7.2,7.5[. For D < 7.2 the 
phase transition is continuous and for T > 7.5 there is no 
phase transition. We show in Fig. [6] (bottom) Tc versus 
D. 

We note that the maximal value of D and the tricritical 
value Dc depend on the film thickness. We can notice 
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FIG. 5: (Color online) Energy (top) and magnetization (bot¬ 
tom) versus MC time t (in unit of 10^) at the transition tem¬ 
perature Tc — 2.694 for D = 7.3 (first-order transition). Note 
that the red and green curves are obtained with ferromagnetic 
and random initial configurations, respectively. See text for 
comments. 




FIG. 6: (Color online) Latent heat AE (top) and Tc (bottom) 
versus D, for L;^ = 4, L = 120. The critical value of D (~ 
7.2) is shown by the vertical arrow. The dotted line between 
D = 7.4 and 7.5 is extrapolated. 


that by looking at the bulk maximal value D = 8 {L^ = 
oo) and Dc = 7.5 as shown in Fig. [TJ The four-layer film 
has Dc = 7.2. So, when Lz goes to infinity Dc goes from 
7.2 to 7.5. 


B. Size effect 


When the system size is infinite, in second-order phase 
transitions the correlation length is infinite at the criti¬ 
cal point. However, in first-order transitions the correla¬ 
tion length is finite at the transition temperature where 
the two phases coexist and the energy is discontinuous. 
In simulations, in spite of the fact that we use periodic 
boundary conditions to mimic large systems, we cannot 
avoid finite-size effects on the results. The nature of the 
transition may not be detected at small system sizes. It 
is therefore very important to measure the size effects in 
numerical simulations. We have shown in Fig. [2] (bot¬ 
tom) the energy versus T for L = 36 and L = 300. The 
size effect is extremely small. The transition remains 
continuous though one observes a change in the slope of 
the curve which is steeper for the larger size. In the first- 
order region, the energy and magnetization are already 
discontinuous even for L as small as 36. 

The film thickness affects, on the other hand, the value 
of the transition temperature Tc as seen in Fig. [71 As 
Lz increases, the transition temperature tends to that of 
the bulk. We have used the least mean-square fit with 
the form: 


T,{L,) = Te(oo) - T ( 7 ) 

Ez 

where A = 2.692 ± 0.165 and Te(oo) = 4.455 ± 0.024. 

The way how Tc increases with increasing thickness 
is characterized by constant A in the above equation of 
Tc{Lz). This constant is different from one material to 
another depending on the coupling between film layers. 
In some materials A is very small, meaning that inter¬ 
layer coupling is very small. This is not the case here in 
spite of the fact that there are only two nearest neighbors 
for each interior atom on the 2 : axis (with only one for 
a surface atom). Knowing how Tc varies with the film 
thickness can help determine the inter-layer coupling. 

At this stage, let us discuss about the criticality of the 
transition in the second-order region. If we compare Eq. 
([7|with the finite-size scaling relation 

T,{L) = Tc(oo) + AL-^!^ (8) 

we see that = 1 which is the 2D Ising universality 
exponent. This is in agreement with Ref. [l^: when the 
film thickness becomes small the critical exponents tend 
to the 2D criticality. 
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FIG. 7: (Color online) Top: Energy versus T for Lz=4: (red), 
8 (green), 12 (blue) and 16 (magenta) (from above), with 
L = 120, and D = 6. Bottom: Tc (blue points) versus Lz for 
D — The continuous red line is the least mean-square fit. 

C. Surface effect 

So far, we have supposed J = K ioi any NN spin pair 
in the film. We investigate now the surface effect due 
to the surface parameter Ks taken to be different from 
K. We write the biquadratic surface and bulk parts as 
follows: 

(9) 

hi 

where ai, = K/ J and ag = Kg/ J denote respectively the 
bulk and surface interactions and ^ denotes the sum 

i'd' 

over NN spin pairs in the surface layer. We take a^, = 1. 
Let us show the magnetization of the first and second lay¬ 
ers in Fig. [8]for several values of ag. As seen, the weaker 
the surface interaction is, the smaller the surface magne¬ 
tization becomes. Only when ag is much larger than 1, 
the surface-layer magnetization becomes larger than the 
second-layer (not shown). For the first-order region, the 
surface and interior magnetizations have discontinuities 
at the transition as expected. 

We show now the average number of spins ±1 and the 
average number of spins zero in each layer versus T in Fig. 
[9] at the first-order transition with D = 7.3. They are 
defined as Mi, 2 (±l) = ^ 6{Si,-1)])/L^ 

where the sum is made for each layer: Mi(±l) [M 2 (±l)] 
corresponds to the surface (second) layer. For spins zero, 
Afi,2(0) = Several remarks are in or¬ 

der: 
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FIG. 8: (Golor online) Layer magnetization of the first (red) 
and second (green) layers versus T. The first three figures 
(from top) are for: as = 0.8, 1, and 1.2, with D = 6. The 
bottom figure represents a first-order case where Os = 1 and 
D = 7.3: red (blue) symbols indicate the first (second)-layer 
magnetization. 


(i) Below the transition temperature, the ordering re¬ 
sults from spins ±1. The number of spins zero increases 
slowly from 0 at T < Tc but becomes dominant for 
T>Tc. 

(ii) At T < Tc the surface has a smaller number of 
spins ±1 than the second layer, namely there is a deficit 
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FIG. 9: (Color online) The normalized numbers of spins ±1 
(He-4) and spins zero (He-3) versus T for the hrst and second 
layers. Red void circles and green circles represent the number 
of He-4 (spins ±1) on the hrst and second layers, while blue 
void squares and magenta squares represent the number of 
He-3 (spins zero) on the hrst and second layers, respectively. 
See text for comments. 

of He-4 at the surface. While, the number of spins zero 
is larger at the surface than in the second layer. 

IV. CONCLUSION 

We have investigated in this paper the BEG model 
used for a thin him of stacked triangular lattices with a 
thickness Lz. There are three important aspects of our 


results for thin hlms: 

(i) the nature of the hrst-order phase transition in a 
region of the phase space is conserved down to a 4-layer 
him, unlike in other systems where bulk hrst-order tran¬ 
sition becomes second-order with small thickness 0, 

(ii) the cross-over from second-order to hrst-order tran¬ 

sition in the bulk is also conserved in thin hlms as shown 
above. The anisotropy of the BEG Hamiltonian affects 
the nature of the phase transition as it has been observed 
in the bulk case of simple cubic lattice in 4-layer 

triangular hlms, for D <1.2 the transition is continuous 
and for 7.2 < D < 7.5 the transition is of hrst order. 
This has been conhrmed with the histogram technique 
where the latent heat can be measured with precision, 

(iii) The surface effect on the layer magnetizations has 
been shown. The surface magnetization is smaller than 
the interior layer if the surface interaction is not so large. 
If we map the BEG model into a mixing of He-3 and He-4, 
then near the surface there is a He-3 enrichment (normal 
liquid) in a him at low temperatures. This point is new 
with respect to the bulk properties where the mixing of 
two liquids is uniform over the system. 
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